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Abstract 

Background: Soybean is one of the most economically important crops in the world. The cotyledon is the nutrient 
storage area in seeds, and it is critical for seed quality and yield. Cotyledon mutants are important for the genetic 
dissection of embryo patterning and seed development. However, the molecular mechanisms underlying soybean 
cotyledon development are largely unexplored. 

Results: In this study, we characterised a soybean curled-cotyledon {ceo) mutant. Compared with wild-type (WT), 
anatomical analysis revealed that the ceo cotyledons at the torpedo stage became more slender and grew outward. 
The entire embryos of ceo mutant resembled the "tail of swallow". In addition, ceo seeds displayed reduced 
germination rate and gibberellic acid (GA 3 ) level, whereas the abscisic acid (ABA) and auxin (IAA) levels were 
increased. RNA-seq identified 1,093 differentially expressed genes (DEGs) between WT and the ceo mutant. 
The KEGG pathway analysis showed many DEGs were mapped to the hormone biosynthesis and signal 
transduction pathways. Consistent with assays of hormones in seeds, the results of RNA-seq indicated auxin 
and ABA biosynthesis and signal transduction in ceo were more active than in WT, while an early step in GA 
biosynthesis was blocked, as well as conversion rate of inactive GAs to bioactive GAs in GA signaling. 
Furthermore, genes participated in other hormone biosynthesis and signalling pathways such as cytokinin 
(CK), ethylene (ET), brassinosteroid (BR), and jasmonate acid (JA) were also affected in the ceo mutant. 

Conclusions: Our data suggest that multiple phytohormone biosynthesis and signal transduction pathways 
are reprogrammed in ceo, and changes in these pathways may partially contribute to the ceo mutant 
phenotype, suggesting the involvement of multiple hormones in the coordination of soybean cotyledon 
development. 
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Background 

Cotyledon is the main nutrient storage area in soybean 
seed, which contains approximately 40% protein and 
20% oil at maturity [1]. Flowering plants are divided into 
monocots and dicots based on the number of cotyle- 
dons. In Arabidopsis, as a paradigm for dicot embryonic 
development, cotyledons are specified from two lateral 
domains at the apical end of the embryo proper at the 
heart stage, showing bilateral symmetry [2]. Several stud- 
ies have reported cotyledon mutants in Arabidopsis [3-7], 
Antirrhinum [8], tomato [9-11], and other plants [12]. 
These mutants have been used to identify a number of 
orthologous gene hierarchies that are involved in cotyle- 
don development [13]. The Arabidopsis cup-shaped coty- 
ledon {CUC) genes (CUC1, CUC2, and CUC3), belonging 
to the plant-specific NAC transcription factor, Picea 
glauca NAC01, and SHOOT MERISTEMLESS (STM) are 
required for shoot apical meristem (SAM) formation 
and cotyledon separation [4,14-16]. STM induces the ex- 
pression of CUC genes, and the induction is specific and 
independent of other meristem regulators [17]. STM nega- 
tively regulates ASYMMETRIC LEAVES1 (AS1), and the 
mutation of AS1 in Arabidopsis disrupts cotyledon devel- 
opment [18]. KNOTTEDAike {KNAT6) is expressed at the 
boundary between the SAM and cotyledon later than 
STM and CUC, suggesting that KNAT6 plays a crucial role 
in SAM maintenance and boundary establishment in em- 
bryos via the STM/CUC pathway [19]. 

Phytohormones, primarily auxin, CK, BR, GA, ABA, 
and JA, orchestrate intrinsic developmental programs. 
Previous studies have greatly advanced the functions of 
individual hormone. During last two decades, extensive 
lists of genes involved in hormone synthesis, catabolism 
and signal transduction have been identified via the ana- 
lysis of mutants [20-23]. Auxin plays a key regulatory 
role in the initiation of cotyledon in the apical margin of 
the globular embryo [24]. In many plant species, includ- 
ing Brassica napus [25], Arabidopsis [26], and Norway 
spruce [27], studies have shown that cotyledon and SAM 
development could be inhibited by various polar auxin 
inhibitors in embryogenesis. PIN auxin efflux regulators 
are the best characterised components of auxin trans- 
port. The polar location of PIN1 is associated with the 
position of auxin accumulation in incipient cotyledons, 
and auxin flows into the primordium interior during 
cotyledon outgrowth [28]. The pinl pid double mutant 
completely lacks cotyledons and bilateral symmetry and 
exhibits increased CUC1, CUC2 and STM expression, 
indicating that directional auxin transport is important 
for the establishment of bilateral symmetry and the pro- 
motion of cotyledon outgrowth [26]. Similarly, the pid 
enp double mutant, which is similar to "laterne", also abol- 
ishes cotyledon development, reflecting failed auxin accu- 
mulation in the apex, and ENP is necessary for cotyledon 



development via the control of PIN1 polarity in the con- 
text of PID [29]. In addition to auxin, some studies have 
shown that CK affects cotyledon development. Exogenous 
application of benzyladenine (the synthetic cytokinin) 
onto hybrid larch somatic embryos reduced the cotyledon 
number [30]. The amp-arabidopsis mutant demonstrated 
polycotyly, indicating an elevated cytokinin level [31]. 
Arabidopsis ESR1 and ESR2, whose expression levels are 
induced by cytokinin, not only control shoot regeneration 
but also play a role in cotyledon development [32,33]. 

In this study, we characterised a soybean mutant named 
ceo, which displayed curled cotyledons. The ceo embryos 
at the torpedo stage resembled the "tail of swallow". 
High performance liquid chromatography (HPLC) analysis 
showed that ceo seeds contained higher IAA and ABA 
levels, but lower GA 3 level. Transcriptome analysis re- 
vealed that multiple phytohormone biosynthesis and sig- 
nal transduction pathways were reprogrammed in ceo. 

Results and discussion 

Embryogenesis in the ceo mutant 

We obtained the ceo mutant following sodium azide 
(NaN 3 ) and 60 Coy ray seed mutagenesis of soybean [34], 
which was initially characterised by outward folding cot- 
yledons (Figure 1A-C, 1M and IN). Compared with WT, 
developing cotyledons of ceo mutant had larger vacuoles 
and more membranous multilamellar structures, as well 
as higher methionine and cysteine content [34]. In 
addition, the germination rates of seeds were compared 
between the ceo mutant and WT. In general, seeds from 
WT germinated on the second day, and the germination 
rate reached -72% on the fifth day. ceo seeds were dor- 
mant until the third day, and the germination rate was 
only -6.7% on fifth day, suggesting that the germination 
rate of ceo seeds was greatly reduced (Figure lO). 

To determine the embryogenesis defects in ceo mutant, 
we traced its phenotype back to the globular stage of em- 
bryogenesis. In higher plants, embryogenesis can be concep- 
tually divided into three overlapping phases: morphogenesis, 
maturation and desiccation [35]. In general, for the wild type 
at the globular stage, an embryo comprises a spherical em- 
bryo proper and a suspensor (Figure ID). The embryo 
proper is spherical during the pro-embryo and globular 
stage and eventually forms cotyledons, shoot meristem 
and hypocotyl region around the apical-basal axis, and this 
structure is connected to the embryo sac via the suspen- 
sor. During the globular-heart transition phase, a dramatic 
change occurs: two of the cotyledon primordia are sym- 
metrically initiated from the lateral domains at the apical 
end of the heart-shaped embryo proper (Figure IE). Sub- 
sequently, the cotyledons and hypocotyl elongate, and the 
embryo enters the torpedo stage in which the cotyledons 
continue to rotate. The morphogenesis phase of embryo- 
genesis ends at the cotyledon stage (Figure IF and 1G). 
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Figure 1 Embryogenesis of the ceo mutant. (A) Cotyledons of WT. (B) Cotyledons of the ceo mutant. (C) Cotyledons after germination for 
5 days. WT (left) and the ceo mutant (right). (D) WT embryo at the globular stage approximately 9 days after fertilisation. (E) WT embryo at the 
heart stage approximately 1 1 days after fertilisation. (F) and (G) WT embryo at the torpedo stage approximately 13 days after fertilisation. (H) ceo 
embryo at the globular stage approximately 1 1 days after fertilisation. (I) ceo embryo at the heart stage approximately 13 days after fertilisation. 
(J) to (L) ceo embryos at the different torpedo stages approximately 15 days after fertilisation. Bar= 100 urn for (D) to (L). (M) WT embryo at the 
cotyledon stage. (N) ceo embryo at the cotyledon stage. (O) Germination rates of WT and ceo. WT, triangles; ceo, circles. **P < 0.01 . 



However, the development of embryos in ceo was 
slower than that in WT (Figure 1). When WT embryos 
were at the heart stage, ceo embryos remained morpho- 
logically at the globular stage (Figure IE and 1H). The 
embryo proper in ceo appeared smaller than that in WT 
(Figure ID and 1H). At the heart stage, cotyledons were 
initiated from the cotyledon primordia. Figure IE and II 
showed that the cotyledon initiation positions were identi- 
cal in both the WT and ceo. However, the WT cotyledons 
grew upwards, while that of the ceo mutant demonstrated 
outward growth. Subsequently, the bifurcation of the coty- 
ledons at the torpedo stage in the ceo mutant was larger, 
and the tips of the cotyledons were sharper (Figure IF to 
1G and 1J to 1L). At the torpedo stage, ceo mutant em- 
bryos resembled "tail of swallow", while the radicle and 
shoot apical meristem were also present and appeared 
normal. All in all, ceo showed abnormal embryogenesis, 
especially at the torpedo stage, when the ceo cotyledons 
became more slender and grew outward. However, SAM 
and the polar axis were normal. Cotyledons and leaves 
share homology in plants [13]. Cotyledon initiation and 
development in ceo were defective, while leaf development 
was correct, suggesting partial homology between leaf and 
cotyledon. 

Increased free IAA and ABA levels, while decreased GA 3 
level in ceo 

The hormone auxin governs a variety of developmental 
processes [36], including embryogenesis [37]. In embryo- 
genesis, auxin controls cell division and specification, 
which are critical for establishing the embryo pattern. 



Most of the pattern formation events in Arabidopsis em- 
bryogenesis, including cotyledon development, depend 
on auxin biosynthesis, transport, and response [37]. IAA 
is the most important natural auxin. Therefore, we ex- 
amined the endogenous free IAA level in WT and the 
ceo mutant seeds at 7 DAF (days after fertilisation) by 
HPLC analysis (Figure 2). Compared with WT, the ceo 
mutant had a higher free IAA level in seeds at 7 DAF. 

In addition, the endogenous GA3 and ABA content in 
WT and ceo seeds at 7 DAF were quantified. The mea- 
surements demonstrated that the ABA content was sig- 
nificantly increased, while the GA3 content was decreased 
in ceo (Figure 2). The increased ABA/GA ratio may lead 
to a low germination rate for the ceo mutant. 
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Figure 2 Quantification of free IAA, GA 3 and ABA in WT and 
the ceo seeds at 7 DAF. **P < 0.01 
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Auxin sensitivity was depressed in ceo 

Twelve auxin responsive genes, including two AUX/IAA 
genes (Glyma08g22190, Glymal0g32340), five GH3 
genes (Glyma02gl3910, Glyma02gl7360, GlymallgOSSlO, 
Glymal2gll890, Glymal2gl 1200) and five SAUR genes 
(Glymal0g06360, Glymal0g06370, Glyma05g36360, Gly- 
mal2g03780, Glyma09g3S460), were selected to investi- 
gate the transcription of early auxin-regulated genes 
in response to auxin treatment by quantitative RT-PCR 
(qRT-PCR). As shown in Figure 3, all of these genes were 
up-regulated in WT at 30 min after treatment with 
0.05 uM 2,4-dichlorophenoxyacetic acid (2,4-D). However, 
in the ceo mutant, only three genes (Glyma08g22190, 
Glymal2gll200, and Glyma09g35460) were weakly 
up-regulated, and the other genes had weak variation or 
even unexpected down-regulation (Figure 3A). When hy- 
pocotyls were treated with 0.2 uM 2,4-D, the twelve genes 
tested showed higher expression levels in WT, and 
Glyma08g22190 and GlymallgOSSlO showed greater 
than five-fold up-regulation (Figure 3B). With the excep- 
tion of Glymal2gll890 and Glyma05g36360, the expres- 
sion trends of the genes in ceo were similar to that 
observed after treatment with 0.05 uM 2,4-D (Figure 3). 
Thus, auxin sensitivity were repressed in the ceo mutant. 



RNA-Seq of WT and the ceo mutant 

The positions of the cotyledon primordia in ceo were gen- 
erally normal, but the abaxial/adaxial patterning of cotyle- 
dons was flawed, which was likely to exist before cotyledon 
initiation. Thus, we focused on the globular stage, and 
the pods at 7 DAF were collected for RNA-Seq using the 
Illumina HiSeq2000 system. A total of 52,870,578 and 
50,494,652 successful reads (average length: -180 bp, ap- 
proximately 4 x soybean genome) were produced for WT 
and the ceo mutant, respectively (Table 1). Notably, more 
than 85% of the reads mapped back to the soybean refer- 
ence genome in Phytozome database [38]. In WT, 80.74% 
of these genes uniquely mapped to a single location, and 
80.27% uniquely mapped to a single location for ceo 
(Table 1). We identified 39,499 highly confident genes in 
WT and 39,527 genes in the ceo mutant. With a False Dis- 
covery Rate (FDR) < 0.001 and |log 2 Ratio| > 1, a total of 
1,093 genes were differentially expressed; 256 DEGs had 
lower expression levels in the ceo mutant, while 837 DEGs 
were activated (Figure 4 and Additional file 1). In this 
study, we primarily focused on the differentially expressed 
genes. The expression levels of the majority of these genes 
were unaffected in the ceo mutant, suggesting that ceo tar- 
geted a limited number of genes. 
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Figure 3 Regulation of early auxin-regulated genes in WT and the ceo muant. (A) Treatment with 5x10 M 2,4-D (B) Treatment with 
2 x 1 0" 4 M 2,4-D. WT, triangles; ceo, circles. 
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Table 1 Summary of RNA-seq data from WT and ceo pods at 7 DAF 





Total reads 


Total base pairs 


Total mapped reads 


Percent mapped 


Unique match 


Percent unique match 


wr 


52,870,578 


4,758,352,020 


45,064,067 


85.23% 


42687787 


80.74% 


ceo 


50,494,652 


4,544,518,680 


42,914,706 


84.99% 


40531216 


80.27% 



To validate the reliability of the expression profiles ob- 
tained using RNA-Seq, we selected fifteen genes with 
higher or lower expression levels for semi-quantitative 
reverse transcription PCR (RT-PCR) analysis (Table 2 
and Additional file 2: Figure S1A and SIB), and three 
(Glyma03g02620, Glymal4g39300, Glymallgll580) of 
the fifteen genes for further qRT-PCR analysis in 7DAF 
pods from WT and the ceo mutant (Figure 5A). For 
all of the genes, the results obtained from the semi- 
quantitative RT-PCR/qRT-PCR analysis were consist- 
ent with data from RNA-Seq. 

In addition, we used semi-quantitative RT-PCR to ana- 
lyse the expression of four selected genes (Glyma03g02620, 
Glymal4g39300, Glymallgll580, Glymal7gl4710) in 
WT and ceo seeds at 7 and 15 DAF (Additional file 2: 
Figure SIC). Using qRT-PCR, we examined the expres- 
sion of three genes {Glyma03g02620, Glymal4g39300, 
Glymallgll580) in the roots, stems, leaves, flowers and 
seeds at the various developmental stages of WT and ceo, 
and the results showed these genes were differentially 
expressed in these tissues (Figure 5B), indicating some of 
the DEGs that we obtained were differentially expressed 
in other tissues. 

Functional classification of differentially expressed genes 

To evaluate the potential functions of the DEGs between 
WT and the ceo mutant, Gene Ontology (GO) categories 
were used to identify key processes for ceo mutant. De- 
tailed GO term annotations of the DEGs were categorised 
according to biological processes, molecular functions 
and cellular components. The GO categories for the set 
of DEGs (Figure 6) revealed that most of the encoded 




Up-regulated Down-regulated 

Figure 4 DEGs between WT and the ceo mutant: the up- and 
down-regulated represent the genes that had an increased and 
decreased expression levels in ceo. 



products were associated with "intracellular organelle", 
"catalytic activity", "response to stimulus" and "cellular 
metabolic process". 

Biological interpretation of the 650 DEGs was further 
completed using KEGG pathway analysis. Many of the 
DEGs were mapped to pathways essential for plant growth 
and development. Overall, 98 pathways were affected in 
the ceo mutant, and some of these pathways were consist- 
ent with biological processes previously revealed in the GO 
analysis. The most frequently represented pathways pri- 
marily involved in secondary metabolism, including plant 
hormone biosynthesis and signal transduction (20.46%), 
such as auxin, CK, GA, ABA, ET, BR and JA (Additional 
file 3). Asakura et al. found that genes expressed in the 
early pod stage of soybean were the most numerous, com- 
pared with 2 and 5 mm seeds, and the number of expressed 
genes gradually decreased with seed development [39]. 
Commonly expressed genes in these stages accounted for 
70-95% of the expressed genes in each stage, and only a 
small number of genes were differentially expressed, which 
mainly included |3-conglycinin, lipid synthesis, lipoxygenase 
and seed-maturation proteins [39]. Therefore, most DEGs 
obtained in this study, including genes that were related to 
hormone biosynthesis and signalling, could still be identi- 
fied in the early seed stage. However, some DEGs may have 
been missed because of tissue-dependent gene expression. 

The results of the transcriptome data showed that many 
genes encoding proteins involved in auxin metabolism and 
signal transduction pathways were differentially expressed 
between ceo and WT. L-tryptophan, as a primary precur- 
sor, is used to generate many indole-containing substances 
in plants, including IAA. a-subunit anthranilate synthase 
(ASA) catalyses the conversion of chorismate to anthra- 
nilate as one of the first steps in the Trp biosynthesis path- 
way. In Arabidopsis, up-regulation of WEAK ETHYLENE 
INSENSITIVE2/ANTHRANILATE SYNTHASE al (WEI2/ 
AS AT) gene through ET results in auxin accumulation in 
the tips of primary roots [40]. Transgenic rice lines ex- 
pressing a-SUBUNIT OF ANTHRANILATE SYNTHASE 
(OASA1D) gene show increased production of tryptophan 
and free IAA [41]. Tryptophan synthase beta subunit (TSB) 
catalyses the last step in tryptophan biosynthesis. In IAA 
biosynthesis, indole-3-acetaldoximem (IAOx)-dependent 
IAA biosynthesis is a specific pathway in plants. In this 
process, two homologous cytochrome P450 enzymes, 
CYP79B2 and CYP79B3, mediate the synthesis of IAOx 
from tryptophan [42], and the double mutation of these 
two genes results in the downregulation of IAA synthesis. 
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Table 2 Fold changes of 15 genes selected from RNA-seq data 



Gene ID 


Gene description 


Fold change 


Regulation 


P-value 


Glyma02g42290 


Transmembrane amino acid transporter protein 


2.02 


Up 


0.000139 


Glyma03g09140 


Transmembrane amino acid transporter protein 


2.69 


Up 


2.58E-07 


Glyma05g23180 


Membrane transport protein 


5.08 


Up 


3.71 E-05 


Glyma08g01430 


WRKY transcription factor 


12.62 


Up 


3.12E-06 


Glyma09g37370 


HLH transcription factor 


2.35 


Down 


2.84E-10 


Glyma16g13570 


bZIP transcription factor 


2.03 


Up 


0.000139 


Glyma18g05720 


NPH3 family 


2.34 


Down 


1.26E-16 


Glyma19g43761 


DUF640 


4.56 


Up 


0.0001 86 


Glyma17g14710 


YABBY protein 


3.05 


Up 


6.86E-24 


Glymal 1g00580 


tryptophan synthase fi-chain 


8.01 


Up 


1.85E-57 


Glyma14g39300 


ubiquitin-protein ligase activity 


371.70 


Up 


3.14E-08 


Glyma03g02620 


LBD transcription factor 


441.39 


Down 


8.88E-05 


Glyma06g02970 


DBB transcription factor 


2.91 


Up 


9.21 E-06 


Glyma06g 17420 


HLH transcription factor 


2.06 


Up 


1 .74E-68 


Glyma16g28310 


Fructose-1 -6-bisphosphatase 


2.16 


Down 


2.39E-43 



In the transcriptome data obtained in this study, one ASA 
gene, one TBS gene and three CYP79 genes were up- 
regulated, suggesting that auxin synthesis was enhanced in 
the ceo mutant (Additional file 4). Mounting evidence sug- 
gests that polar auxin transport (PAT) controls important 



growth and developmental processes in higher plants. In 
Arabidopsis, studies have shown that the AUXIN1/LIKE- 
AUX1 (AUX/LAX) family of auxin transporters are major 
influx carriers, whereas the PIN-FORMED (PIN) family of 
auxin proteins are major efflux carriers. The AUX/LAX 
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Figure 5 Confirmation of RNA-seq results using qRT-PCR. Three DEGs (Glymal 4g39300, Glymal 1g00580, Glyma03g02620) were selected for 
confirmation. (A) Expression levels of three DEGs in pods at 7 DAF. (B) Expression levels of three DEGs in roots, stems, leaves, flowers and seeds 
at various stages (7, 10, 15, 30 and 40 DAF) of WT and ceo. 
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family is represented by four highly conserved genes, 
namely AUX1, LAX1, LAX2, and LAX3. LAX2 regulates 
vascular patterning in cotyledons. Transcriptome data 
showed that three LAX2 genes and one PIN gene were 
up-regulated in the ceo mutant, indicating that auxin polar 
transport was also influenced. Auxin-responsive genes in- 
clude GH3 genes, small auxin up RNA (SAUR) genes and 
AUX/IAA genes [43]. In our transcriptome data, one 
AUX/IAA and five SAUR genes were down-regulated, and 
four SAUR genes were up-regulated in ceo (Figure 7). All 
of the GH3 genes were up-regulated, which were likely 
due to the down-regulation of the AUX/IAA genes as 
Aux/IAA proteins repress the expression of GH3 genes in 
soybean [44]. Taken together, these data indicated that 
auxin biosynthesis and signal transduction in ceo were 
more active than in WT. 

The role of GA in embryo morphogenesis has not yet 
been determined. In general, GA is required for seed de- 
velopment, which is determined through the analysis of 
GA-deficient mutants. In pea {Pisum sativum L.), the GA- 
biosynthesis mutant Ih has provided the best evidence for 
the physiological role of GA in seed development. The LH 
protein, which encodes e«f-kaurene oxidase (KO), is 
required for the three-step oxidation of e«i-kaurene to 
e«f-kaurenoic acid as an early step in the GA biosynthesis 
pathway altering seed development [45]. KS (e«t-kaurene 
synthase) and KO are key enzymes in the GA biosynthetic 
pathway, and both are encoded by a single gene in soy- 
bean. GA 2 ox (GA 2 -oxidase) plays an important role in the 
last step of GA biosynthesis and catalyses the inactivation 
of GAs by conversion to bioactive GAs, thereby affecting 
the amount of bioactive GAs. Overexpression of GA 2 ox 
genes results in reduced GA level in Arabidopsis [46], rice 
[47], and poplar [48]. DELLA-domain proteins are tran- 
scriptional regulators that repress GA responses, and 



these proteins are rapidly degraded in response to GA In 
this study, one KO gene and one KS gene were downregu- 
lated in the ceo mutant, suggesting that an early step in 
GA biosynthesis was blocked (Additional file 5). In 
addition, four GA^px genes and three DELLA genes were 
up-regulated, suggesting that the rate of conversion of 
inactive GAs to bioactive GAs was limited for the ceo 
mutant (Figure 7 and Additional file 5). 

ABA not only plays an important role in the stress re- 
sponses and plant tolerance, but also controls seed develop- 
ment and germination. The enzyme 9-cis-epoxycarotenoid 
dioxygenase (NCED) is involved in a rate-limiting step for 
ABA biosynthesis. Recently, Martfnez-Andujar et al. re- 
ported that the NCED6 induction enhanced seed dor- 
mancy in Arabidopsis [49]. Here, we showed that NCED 
genes were up-regulated in the ceo mutant, suggesting that 
ABA synthesis was promoted. However, another gene 
encoding abscisic acid 8 -hydroxylase, which catalyses 
a key step in ABA catabolism, was also up-regulated 
(Additional file 6). Snfl-related protein kinases (SnRK2s) 
are key positive regulators in ABA signal transduction 
pathway. In the absence of ABA, protein phosphatases 
type 2C (PP2Cs) physically interact with SnRK2s and 
negatively regulate ABA responses. However, interactions 
between PP2Cs and SnRK2s are disrupted when an ABA 
molecule binds to PYRABACTIN RESISTANCE1 (PYR1)/ 
PYR1-LIKE (PYL)/REGULATORY COMPONENTS OF 
ABA RECEPTORS (RCAR) receptors, leading to struc- 
tural changes in these receptors. When SnRK2s are 
released from PP2C inhibition, their downstream tar- 
gets, including ABA-Responsive Element Binding Factors 
(ABFs), are activated [50]. Fujii et al. constructed an 
Arabidopsis line carrying mutations in all 10 members of 
the SnRK2 family. The decuple mutant snrk2. 1/2/3/4/5/6/ 
7/8/9/10 (srk2g/d/e/a/h/e/f/c/j/b) was defective in gene 
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regulation and ABA accumulation [51]. In our data, the 
PP2Cs, SnRK2 and ABF genes were up-regulated in the 
ceo mutant (Figure 7). In summary, ABA biosynthesis, ca- 
tabolism and signalling were enhanced in ceo. 

In addition to auxin, GA and ABA, the pathways of 
other hormone biosynthesis and signal transduction 
(such as CK, ET, BR, and JA) were also influenced in the 
ceo mutant. CK and ET, as major classes of plant hor- 
mone, are involved in various aspects of plant develop- 
ment, including organ development [52,53]. In ceo, we 
found two genes encoding adenosine-phosphate isopente- 
nyltransferase (IPT), a rate-limiting enzyme in CK biosyn- 
thesis, expressed at significantly higher levels [54,55] 
(Additional file 7). In CK signalling, type-B ARABIDOPSIS 
RESPONSE REGULATORS (ARRs) are transcriptional ac- 
tivators that regulate cytokinin targets, including type-A 
ARRs [56]. We observed the up-regulation of two type-B 
ARR genes (Figure 7). Aminocyclopropane-carboxylate 
oxidase (ACO), an "ethylene-forming enzyme" (EFE), 
catalyses the final step in ethylene biosynthesis, con- 
verting 1-amino-cyclopropane-l-carboxylic acid (ACC) 
into ethylene, C02, and cyanide. Four ACO genes were 
up-regulated in ceo (Additional file 8). Antisense con- 
structs of ACO genes can notably reduce ET production 
[57]. ET signalling was also influenced in the ceo mutant, 
and the expression levels of ethylene-related genes en- 
coding serine/threonine-protein kinase (CTR1), mitogen- 
activated protein kinase 6 (MPK6), and ethylene-responsive 
transcription factor (ERF1/2) were altered in ceo (Figure 7 
and Additional file 8). BRs regulate embryonic and post- 
embryonic development, adult homeostasis and expression 
of numerous genes [58]. Plant steroid hormones are recog- 
nised by brassinosteroid-insensitive 1 (BRI1), whose kinase 
activity is activated by exogenous BR application [59]. We 
also found that the expression levels of thirteen BR1 and 
four BRI1 Associated receptor Kinase 1 (BAK1) genes in ceo 
were higher than those in WT (Figure 7). J A regulates plant 
growth and stress responses, biosynthetic pathway of which 
starts with oc-linolenic acid [60]. The DEGs in JA biosyn- 
thesis and signalling pathway included those that encode 
phospholipase Al, 12-oxo-phytodienoic acid reductase, 
lipoxygenase, JASMONATE-ZIM (JAZ) and MYC2 pro- 
teins [61] ( Additional file 9 and Figure 7). 

Auxin, BR, CK and GA are generally considered major 
developmental growth regulators, whereas ABA, ET and 
JA are involved in stress responses, all of these hor- 
mones converge on auxin [62,63]. Auxin is not only a 
hormone but also a morphogen. Auxin accumulation is 
controlled by auxin biosynthesis and metabolism, the 
conjugation/deconjugation of active auxin to/from its in- 
active conjugated form and auxin transport [64-66]. 
Auxin maxima and activity have been implicated in em- 
bryo development and patterning [37]. In ceo, the free 
IAA level and auxin polar transport were affected. Polar 



auxin transport promotes the formation of local auxin 
maxima and gradients within tissues that play important 
roles in auxin action, such as the patterns of cell division 
and differentiation in the root meristem [67] and the es- 
tablishment of cotyledons [26,29]. Plants show direc- 
tional growth in response to environmental stimuli such 
as light or gravity (tropisms). It was found that auxin 
played an important role in tropisms [68] and unequal 
auxin distribution across organs may lead their bending 
[69]. The cotyledons of ceo showed outward bending, 
which allowed us to hypothesise that PAT might cause 
an auxin imbalance between the abaxial and adaxial side 
in cotyledons of the ceo mutant, resulting in the adaxial 
side developing faster than the abaxial side. However, in 
soybean embryogenesis, the endosperm surrounds the 
embryo and cotyledons, the curled cotyledons observed in 
the ceo mutant might reflect the fetter of the endosperm. 
Many signals can regulate auxin polar transport including 
stress, protein phosphorylation, reactive oxygen, flavonols 
and other hormones [70-73], which among these signals 
affect(s) PAT in ceo require further investigation. In ceo, 
only the shape of cotyledons not including leaves was ab- 
normal, indicating that signal(s) or mutated genes regulat- 
ing PAT were specific for cotyledon initiation and growth. 

GA and ABA have mainly been implicated in later 
stages of embryo development, acting antagonistically in 
the regulation of seed germination [74]. Compared with 
WT, ceo showed higher IAA and ABA levels, whereas a 
lower GA 3 level. Several reports indicate that auxin and 
GA overlap in the regulation of multiple aspects of plant 
development, and there is positive cross-talk interaction 
between them. Willige et al. have reported that a mutant 
deficient in GA biosynthesis and signalling showed a re- 
duced PIN protein level, which could be recovered by 
exogenous GA application. Moreover, GA could interact 
with PIN1 and promote cotyledon differentiation [75]. 
Auxin can also induce GA biosynthesis. Root cells with 
lower auxin level have reduced GA biosynthesis [76,77]. 
Interestingly, our results demonstrated a lower GA 3 but 
higher IAA level in ceo seeds. But beyond that, ceo seeds 
had higher ABA content, which might correlate with 
lower GA 3 level and may be the reason for lower ger- 
mination rate in ceo. Consistently, ABA negatively regu- 
lates bioactive GA level during seed germination. 

Conclusions 

Multiple phytohormone biosynthesis and signalling path- 
ways are reprogrammed in the ceo mutant, which is con- 
sistent with our HPLC analysis. Our results indicate that 
curled-cotyledons and reduced seed germination rate may 
result from the biosynthesis and signalling disorder of 
multiple phytohormones. Accumulating evidence sup- 
ports the idea that hormones not only act in a linear path- 
way but also function in a complex network of feedback 
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regulation and cross talk within or between different 
hormone pathways. At present, the studies of molecular 
mechanisms underlying hormonal cross talk mainly focus 
on postembryonic development, including the roots, 
shoots and leaves. However, adequate models of the hor- 
monal network involved in cotyledon development are 
not well established in plants. Analysis of ceo may provide 
a promising frame for further studies of hormone interac- 
tions and offer the potential discovery of new genes con- 
trolling cotyledon development in plants. 

Methods 

Plant materials 

Soybean seeds from wild type cultivars Nannong 94-16 
and its cotyledon mutant ceo were provided by Soybean 
Research Institute, Nanjing Agricultural University, China. 
The seeds were cultivated in field under natural condi- 
tions in Nanjing, China. Total RNA was extracted from 
the WT and ceo tissues, including the roots, stems, leaves, 
flowers and seeds, which were frozen in liquid nitrogen 
and stored at -80°C until further analysis. 

Plant paraffin section 

For paraffin sectioning, wild type and mutant seeds at 
various developmental stages were fixed in Carnoy's fluid 
(ethanol/glacial acetic acid, 3:1) at room temperature 
followed by dehydration in a graded ethanol series, 
staining with 1% haematoxylin solution, and clearance 
with xylene. The seeds were embedded in paraffin and 
sectioned at 10 to 15 um (Leica, RM2135). 

Semi-quantitative RT-PCR and quantitative RT-PCR 
(qRT-PCR) 

To analyse the gene expression in soybean, total RNA 
was extracted with the Plant RNA Extract Kit (TianGen, 
Beijing, China) according to the manufacturer's instruc- 
tions. PCR amplification was performed with PrimeScript 
1 st Strand cDNA Synthesis Kit (Takara). qRT-PCR was 
performed as described previously [78]. The Gmtubulin 
(GenBank no.AY907703) and Gmactin (GenBank no. 
XM 003531354) genes were used as reference genes 
for qRT-PCR and semi- quantitative RT-PCR, respect- 
ively. The primers used in this study were provided 
in Additional file 10. 

Extraction of endogenous free IAA, GA 3 , ABA and HPLC 
analysis 

A total of 0.5 g fresh seeds at 7DAF from at least 10 in- 
dividual lines of WT and ceo were ground in liquid ni- 
trogen, phytohormones were extracted in 10 milliliter 
cold 80% (v/v) methanol and 1% (w/v) polyvinylpyrroli- 
done was added as an antioxidant. After 24 h of incuba- 
tion at 4°C, the samples were centrifuged for 15 min at 
12,000 g at 4°C. The supernatant was retained, and the 



pellet was re-extracted with 10 ml of cold 80% (v/v) 
methanol and recentrifuged as indicated above. The su- 
pernatants were pooled. After evaporation under a N2 
gas flow, the pH of the remaining water phase was ad- 
justed to 3.0. IAA, GA 3 and ABA purification was per- 
formed using acetyl acetate. The acetyl acetate was also 
evaporated under an N2 stream. The resulting dried pre- 
cipitate was collected in 0.5 ml of 100% methanol, filtered 
through 0.45 um filter membrane and submitted to HPLC 
(Shimadzu, Kyoto, Japan) analysis. HPLC analysis was per- 
formed using an Alltima C18 column (4.6 mm x 250 mm). 
The IAA and ABA concentration was determined at 
254 nm, and the GA 3 concentration was determined at 
210 nm. Hormone measurements were conducted in trip- 
licate and subjected to statistical analysis. 

IAA sensibility assay 

WT and ceo seeds were germinated in moistened ver- 
miculite at 30°C. After 72 h, the seedlings were sprayed 
with 2.5 x 10" 3 M 2,4-D solution (pH 6). The intact 
hypocotyls were then segregated from the seedlings and 
incubated in 10 ml of incubation medium in the absence 
or presence of 5 x 10 s M and 2 x 10~ 4 M 2,4-D (pH 6) 
at 30°C with continuous shaking. At 0.5 and 1 h, the 1.2- 
cm sections obtained from the elongating region of the 
hypocotyls (0.5-1.7 cm below the cotyledons) were ex- 
cised. The hypocotyl sections were rinsed in deionised 
distilled H 2 0 and immediately frozen in liquid N 2 . The 
sections were stored at -70°C until further use for RNA 
extraction. 

cDNA library construction for RNA-Seq and sequencing 

Beads with oligo (dT) were used to isolate poly (A) mRNA 
after total RNA was collected. Fragmentation buffer was 
added to digest the mRNA to generate short fragments. 
These short fragments were used as templates, and ran- 
dom hexamer- primers were used to synthesise first-strand 
cDNA. Second-strand cDNA was synthesised in a reaction 
containing buffer, dNTPs, RNase H and DNA polymerase 
I. Short fragments were purified using the QIAquick PCR 
extraction kit (QIAGEN) and eluted in EB buffer for end 
preparation and poly (A) addition. Subsequently, the short 
fragments were ligated with sequencing adaptors. For 
PCR amplification, we selected suitable fragments to use 
as templates based on the results of agarose gel electro- 
phoresis. The libraries were sequenced using the Illumina 
HiSeq™ 2000 system. 

Raw read filtering and mapping to the reference genome 
and genes sequences 

Dirty raw reads were discarded via the following four 
steps: 1. reads with adaptors were removed; 2. reads with 
unknown nucleotides larger than 5% were removed; 3. 
reads with low quality (more than half of the base quality 
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less than 5) were removed; 4. clean reads were obtained. 
Clean reads were mapped to the reference genome and 
genes sequences respectively using SOAPaligner/soap2. 
Mismatches with no more than m (m default is 5) bases 
were used in the alignment. All subsequent analyses were 
based on clean reads. 

Gene expression annotation 

Gene coverage is the percentage of a gene covered by 
reads. This value equals the ratio of the number of bases 
in a gene covered by unique mapping reads to the num- 
ber of total bases in that gene. The UniGene expression 
was calculated using the RPKM method (Reads Per kb 
per Million reads). The RPKM method eliminated the 
influence of different gene lengths and sequencing dis- 
crepancies on the gene expression calculations. There- 
fore, the calculated gene expression could be used to 
directly compare differences in gene expression between 
the samples. If there was more than one transcript for a 
gene, the longest transcript was used to calculate the ex- 
pression level and coverage. 

Differentially expressed gene analysis 

Using "The significance of digital gene expression pro- 
files" [79], we identified differentially expressed genes 
between WT and the ceo mutant based on the following 
criteria: (FDR) < 0.001 and |log 2 Ratio| > 1. 

GO enrichment and KEGG pathway analysis 

The differentially expressed genes were subjected to GO 
enrichment analysis. The analysis first mapped all DEGs to 
GO terms in the database (http://www.geneontology.org/), 
calculating gene numbers for every term, followed by the 
ultra-geometric test to identify significantly enriched GO 
terms in DEGs compared with the genome background. 
KEGG Pathway (http://www.genome.jp/kegg/pathway.html) 
enrichment analysis was used to identify significantly 
enriched metabolic or signal transduction pathways in 
DEGs compared with the whole genome background. 

Availability of supporting data 

The RNA-seq data of this article have been deposited in 
the Gene Expression Omnibus (GEO) Database, with the 
following access number: GSE58354 (http://www.ncbi.nlm. 
nih.gov/geo/query/acc.cgi?acc=GSE58354). Other support- 
ing data are included as additional files. 

Additional files 



Additional file 1: 1093 differentially expressed genes between WT 
and ceo. 

Additional file 2: Semi-quantitative RT-PCR analysis of 15 selected 
genes. (A) and (B): semi-quantitative RT-PCR analysis of 15 selected 
genes in WT and ceo pods at 7 DAF. (C): semi-quantitative RT-PCR 



analysis of 4 of the 15 selected genes in WT and ceo seeds at 7 and 15 
DAF. 1 : WT seeds at 7 DAF; 2: ceo seeds at 7 DAF; 3: WT seeds at 1 5 DAF; 
4: ceo seeds at 1 5 DAF. 

Additional file 3: Differentially expressed genes in hormone 
biosynthesis and signal transduction pathways. 

Additional file 4: Transcriptional changes in tryptophan 
metabolism. Differentially expressed genes from RNA-Seq were mapped 
to the KEGG pathway database. The up-regulated genes in ceo are boxed 
in green. 

Additional file 5: Transcriptional changes in zeatin biosynthesis. 

The up-regulated and down-regulated genes in ceo are boxed in green 
and red, respectively. 

Additional file 6: Transcriptional changes in diterpenoid 
biosynthesis. The up-regulated and down-regulated genes in ceo are 
boxed in green and red, respectively. 

Additional file 7: Transcriptional changes in carotenoid 
biosynthesis pathway. The up-regulated and down-regulated genes in 
ceo are boxed in green and red, respectively. 

Additional file 8: Transcriptional changes in cysteine and methionine 
metabolism pathway. The up-regulated and down-regulated genes in 
ceo are boxed in green and red, respectively. 

Additional file 9: Transcriptional changes in a-linolenic acid 
metabolism pathway. The up-regulated and down-regulated genes 
in ceo are boxed in green and red, respectively. 

Additional file 10: The list of primers used in this article. 
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